"""
基于19.8年的卫星测高轨迹点数据进行最小二乘法调和分析，得到调和常数。
生成的gridtiff格式的调和常数文件储存方式为：
   波段0: K2_H      波段1: K2_G
   波段2: N2_H      波段3: N2_G
   波段4: S2_H      波段5: S2_G
   波段6: M2_H      波段7: M2_G
   波段8: Q1_H      波段9: Q1_G
   波段10: O1_H     波段11: O1_G
   波段12: K1_H     波段13: K1_G
   波段14: P1_H     波段15: P1_G
"""

import time
import datetime
import math
from collections import deque
from multiprocessing import Pool
import os
import sys
import netCDF4 as nc
import numpy as np
import pandas as pd

root_path = r"__file__"                                                                     # 添加pydggs路径
if root_path not in sys.path:
    sys.path.insert(0, root_path)

from PyDGGS.PyDGGS import *
from common import NAN, MAX_LONGTITUDE, ONE_HOUR_TO_SECOND, HOURS_PER_DAY, level, grid_type, element_type, σ, tide_start_time, cal_f_u, cal_v    # 导入常量

NC_DATA_START_DATETIME = datetime.datetime(1950, 1, 1)                                      # 测高数据的时间单位为自1950年1月1日起n天
TRACK_CYCLE = datetime.timedelta(days=9.9156)                                               # 卫星轨迹周期

input_dir = os.path.join(os.path.dirname(r'__file__'))                                      # 数据位置
input_elev_path = os.path.join(input_dir, 'input', 'SeaLevel')                              # 输入卫星测高数据
ocean_grid_path = os.path.join(input_dir, "output", "Ocean_Grid","Ocean_Grid.json")         # 海洋格网文件
temp_output_path = os.path.join(input_dir, 'temp', 'tide')                                  # 临时工作目录，存储从原始nc文件中提取数据到以格网编码为索引的CSV数据文件。
hg_output_path = os.path.join(input_dir, 'output', 'hg', 'hg.gridjson')                     # 存储提取出来的调和常数


def createIndexCsv(level, grid_type, element_type, netcdf_file, index_file):
    """
    创建索引文件

    :param level:        离散格网等级
    :param grid_type:    离散格网类型
    :param element_type: 离散格网元素类型
    :param netcdf_file:  **NetCDF**所在路径
    :param index_file:   索引文件**CSV**的输出路径
    """
    py_kernel = CreateKernelObj(grid_type)

    # 打开NetCDF文件
    with open(netcdf_file, mode='rb') as f:
        with nc.Dataset(netcdf_file, mode='r', memory=f.read()) as ncfile:

            # 读取变量
            ocean_tide_nc = ncfile["ocean_tide"][:]
            ocean_tide_nc[np.isnan(ocean_tide_nc)] = NAN
            time_nc = ncfile["time"][:]
            lon_nc = ncfile["longitude"][:]
            lat_nc = ncfile["latitude"][:]

    lon_max = np.max(lon_nc)

    data_list = []
    for i in range(time_nc.shape[0]):

        time = time_nc[i]
        lat = lat_nc[i]
        lon = lon_nc[i]
        if lon_max > MAX_LONGTITUDE and lon > MAX_LONGTITUDE:
            lon = lon - 360
        ocean_tide = ocean_tide_nc[i]

        code = DGGSElementCode()
        code.elementType = element_type
        code.level = level
        coord = GeoCoord()
        coord.lat = lat
        coord.lon = lon
        py_kernel.encode(coord, code)

        code_id = code2Id(code, grid_type)

        new_data = {
            'id': code_id,
            'datatime': time,
            'ocean_tide': ocean_tide
        }

        data_list.append(new_data)

    DestroyKernelObj(py_kernel)

    df = pd.DataFrame(data_list)

    df = df.sort_values(by='id')

    df.to_csv(index_file, index=False)
    print(index_file, "生成结束")


def createIndexPath(level, grid_type, element_type, netcdf_path, index_path):
    """
    创建索引文件夹，多进程运行

    :param level:        离散格网等级
    :param grid_type:    离散格网类型
    :param element_type: 离散格网元素类型
    :param netcdf_path:  **NetCDF**所在目录
    :param index_path:   索引文件**CSV**的输出目录
    """
    items = os.listdir(netcdf_path)

    netcdf_files = [
        os.path.splitext(item)[0]
        for item in items
        if os.path.isfile(os.path.join(netcdf_path, item))
    ]

    if not os.path.exists(index_path):
        os.makedirs(index_path)

    tasks = []
    for netcdf_file in netcdf_files:
        netcdf_file_path = os.path.join(netcdf_path, netcdf_file + ".nc")
        index_file_path = os.path.join(index_path, netcdf_file + ".csv")
        tasks.append((level, grid_type, element_type, netcdf_file_path, index_file_path))

    p = Pool(4)
    for task in tasks:
        p.apply_async(createIndexCsv, task)

    p.close()
    p.join()


def calLse(wave_names, tide_list, delta=None):
    """
    最小二乘法来源于陈宗镛的潮汐学第六章的方法
    最小二乘法：least squares estimate (lse)
    原理：实测潮位real_tide 与 计算所得潮位predic_tide 的误差平方和最小
    特点：求系数ai时与系数bi无关
    时间范围：时间间隔1小时，时间长度369，即T = 369 * 24
    所求系数[a0,a1,a2,....a8] a0为截距  [b1,b2,.....b8]

    :param wave_names:   分潮名称List
    :param tide_list:    按时间排序的潮位 [(潮位1, 时间1), (潮位2, 时间2), ...]
    :param delta:        相邻时刻的差，在其时间邻域内的相邻两个潮位被视为同一个点的两个时刻的潮位
    :return:             所求系数[a0,a1,a2,....a8] a0为截距  [b1,b2,.....b8]
    """
    if not tide_list or len(tide_list) == 0:
        return None
    if not delta:
        delta = 1 * TRACK_CYCLE
    min_delta = delta - delta / 16  # 在这个时间上下限内取值，在这里经过测试取16
    max_delta = delta + delta / 16

    tide_elevations = [tide_list[0][0]]
    tide_elevations_length = len(tide_elevations)
    for tide_index in range(tide_elevations_length):
        if isinstance(tide_elevations[tide_index], str):
            tide_elevations[tide_index] = float(tide_elevations[tide_index])

    start_datetime = tide_list[0][1]
    missing_value_indexs = []
    for elevation, cur_datetime in tide_list[1:]:
        if cur_datetime >= start_datetime + max_delta:
            if len(tide_elevations) == 1:
                return None

            while cur_datetime >= start_datetime + max_delta:
                missing_value_indexs.append(len(tide_elevations))
                tide_elevations.append(tide_elevations[-1])
                start_datetime += delta

            tide_elevations.append(elevation)
            start_datetime = cur_datetime

        if cur_datetime >= start_datetime + min_delta:
            tide_elevations.append(elevation)
            start_datetime = cur_datetime

    a = None
    b = None

    wave_count = len(wave_names)
    last_waves_HG = None
    for iteration in range(6):  # 进行多次循环，迭代进行调和分析，以弥补轨迹点缺失值带来的误差，经测试，6次效果最佳
        T = len(tide_elevations)
        if T < 100:  # 100  #每个格元至少有100个数据，才进行调和分析，为了提取的精度。
            return None

        α = np.zeros((wave_count + 1, wave_count + 1))
        for i in range(wave_count + 1):
            for j in range(wave_count + 1):
                if i == j:
                    if i == 0:
                        α[i][j] = 2
                    else:
                        α[i][j] = 1 + math.sin(math.radians(
                            σ[wave_names[i - 1]] * T)) / (σ[wave_names[i - 1]] * T)
                else:
                    if i == 0:
                        α[i][j] = 2 * math.sin(math.radians(
                            σ[wave_names[j - 1]] * T / 2)) / (σ[wave_names[j - 1]] * T / 2)
                    elif j == 0:
                        α[i][j] = 2 * math.sin(math.radians(
                            σ[wave_names[i - 1]] * T / 2)) / (σ[wave_names[i - 1]] * T / 2)
                    else:
                        temp_1 = math.sin(math.radians((
                                                               σ[wave_names[i - 1]] - σ[wave_names[j - 1]]
                                                       ) * T / 2)) / (
                                             (σ[wave_names[i - 1]] - σ[wave_names[j - 1]]) * T / 2)
                        temp_2 = math.sin(math.radians(
                            (σ[wave_names[i - 1]] + σ[wave_names[j - 1]]) * T / 2
                        )) / ((σ[wave_names[i - 1]] + σ[wave_names[j - 1]]) * T / 2)
                        α[i][j] = temp_1 + temp_2

        α_y = np.zeros((wave_count + 1, 1))
        for i in range(wave_count + 1):
            if i == 0:
                sum_0 = 0
                for k in range(T):
                    sum_0 += float(tide_elevations[k])
                α_y[i, 0] = 2 / T * sum_0
            else:
                value_sum = 0
                t = - T / 2 * delta.total_seconds() / ONE_HOUR_TO_SECOND
                for k in range(T):
                    value_sum += tide_elevations[k] * math.cos(math.radians(
                        σ[wave_names[i - 1]] * t
                    ))
                    t = t + delta.total_seconds() / ONE_HOUR_TO_SECOND
                α_y[i, 0] = 2 / T * value_sum

        a = np.linalg.solve(α, α_y)

        β = np.zeros((wave_count, wave_count))
        for i in range(wave_count):
            for j in range(wave_count):
                if i == j:
                    β[i][j] = 1 - math.sin(math.radians(
                        σ[wave_names[i]] * T)) / (σ[wave_names[i]] * T)
                else:
                    temp_1 = math.sin(math.radians((
                                                           σ[wave_names[i]] - σ[wave_names[j]]) * T / 2
                                                   )) / ((σ[wave_names[i]] - σ[wave_names[j]]) * T / 2)
                    temp_2 = math.sin(math.radians((
                                                           σ[wave_names[i]] + σ[wave_names[j]]) * T / 2
                                                   )) / ((σ[wave_names[i]] + σ[wave_names[j]]) * T / 2)

                    β[i][j] = temp_1 - temp_2

        β_y = np.zeros((wave_count, 1))
        for i in range(wave_count):
            value_sum = 0
            t = - T / 2 * delta.total_seconds() / ONE_HOUR_TO_SECOND
            for k in range(T):
                value_sum += tide_elevations[k] * math.sin(math.radians(σ[wave_names[i]] * t))
                t = t + delta.total_seconds() / ONE_HOUR_TO_SECOND
            β_y[i, 0] = 2 / T * value_sum

        b = np.linalg.solve(β, β_y)

        min_hours = (tide_list[0][1] - tide_start_time).total_seconds() / ONE_HOUR_TO_SECOND
        max_hours = (tide_list[-1][1] - tide_start_time).total_seconds() / ONE_HOUR_TO_SECOND
        (a0, waves_HG) = calHG((a, b), wave_names, (min_hours + max_hours) / 2)

        if last_waves_HG:
            for wave_name in waves_HG:
                if abs(last_waves_HG[wave_name]['H'] - waves_HG[wave_name]['H']) < 1e-8 and abs(
                        last_waves_HG[wave_name]['G'] - waves_HG[wave_name]['G']) < 1e-8:
                    break
        last_waves_HG = waves_HG

        for index in missing_value_indexs:
            elevation_datetime = tide_list[0][1] + index * delta
            elevation = calElevation(waves_HG, wave_names, elevation_datetime)

            tide_elevations[index] = a0 + elevation

    return (a, b)


def calHG(lse_result, wave_names, hours):
    """
    计算调和常数

    :param lse_result:  最小二乘法的结果, 包括：最小二乘法计算a0, a1 ~ a8, b1 ~ b8
    :param wave_names:  分潮名称
    :param hours:       潮位时间序列的中间时刻
    :return:            不同分潮的hg
    """
    wave_V = cal_v(hours)
    wave_fu = cal_f_u(hours)

    wave_count = len(wave_names)
    wave_hg = {}
    a0 = float(lse_result[0][0])
    for wave_index in range(wave_count):
        wave_name = wave_names[wave_index]

        a = lse_result[0][wave_index + 1]
        b = lse_result[1][wave_index]

        R = math.sqrt(a * a + b * b)
        Θ = math.degrees(math.atan(b / a))

        if b >= 0 >= a:
            Θ = 180 - abs(Θ)
        elif b <= 0 and a <= 0:
            Θ = abs(Θ) + 180
        elif a <= 0 <= b:
            Θ = 360 - abs(Θ)

        wave_fcosu = wave_fu[wave_name]['fcosu']
        wave_fsinu = wave_fu[wave_name]['fsinu']
        wave_f = math.sqrt(wave_fcosu * wave_fcosu + wave_fsinu * wave_fsinu)
        angle = math.degrees(math.atan(wave_fsinu / wave_fcosu))

        if wave_fsinu >= 0 and wave_fcosu <= 0:
            wave_u = 180 - abs(angle)
        elif wave_fsinu <= 0 and wave_fcosu <= 0:
            wave_u = abs(angle) + 180
        elif wave_fsinu <= 0 and wave_fcosu >= 0:
            wave_u = 360 - abs(angle)
        else:
            wave_u = angle

        H = R / wave_f
        G = (Θ + wave_u + wave_V[wave_name]) % 360

        wave_hg[wave_name] = {'H': H, 'G': G}

    return (a0, wave_hg)

def id2CodeToObj(code_id, grid_type, level, element_type):
    """
    将ID转换为DGGSElementCode对象

    :param code_id:       元素ID
    :param grid_type:     网格类型
    :param level:         网格等级
    :param element_type:  元素类型
    :return: code         DGGSElementCode对象
    """
    code = DGGSElementCode()
    code.level = level
    code.elementType = element_type
    id2Code(code_id, code, grid_type)
    return code


def saveHgFile(level, grid_type, element_type,
                netcdf_path, index_path, ocean_grid_file, out_file,
                is_interpolate=False, interpolation_count=0):
    """
    获取gridtiff格式的hg文件

    :param level:               离散格网等级
    :param grid_type:           离散格网类型
    :param element_type:        离散格网元素类型
    :param netcdf_path:         输入的卫星轨迹点数据**NetCDF**所在目录
    :param index_path:          索引文件**CSV**的输出目录
    :param ocean_grid_file:     海洋格网路径
    :param out_file:            **GridTIFF**格式的振幅和迟角文件保存路径
    :param is_interpolate:      是否对振幅和迟角插值
    :param interpolation_count: 插值次数
    """
    # 潮位、时刻
    wave_names = ['K2', 'N2', 'S2', 'M2', 'Q1', 'O1', 'K1', 'P1']  # 八大分潮

    # 生成索引表
    createIndexPath(level, grid_type, element_type, netcdf_path, index_path)

    # 创建一个空的 DataFrame 作为索引表
    index_df = pd.DataFrame()
    # 读取索引表并按编码、时间排序
    all_files = [os.path.join(index_path, f) for f in os.listdir(index_path) if f.endswith('.csv')]
    number = len(all_files)
    finish_number = 0
    for file in all_files:
        df = pd.read_csv(file)
        df.set_index(['id', 'datatime'], inplace=True)  # 设置多级索引

        if finish_number == 0:
            index_df = df
        else:
            index_df = pd.concat([index_df, df])
        finish_number += 1

        print(finish_number, number)
    index_df = index_df.sort_values(by='id')

    # 获取所有不重复的 id
    unique_ids = index_df.index.get_level_values('id').unique()
    # 逐格元计算其HG并保存
    result_num = 0

    kernel = CreateKernelObj(grid_type)

    # 加载海洋格网
    ocean_layer = DGGSFieldLayerInt()
    ocean_layer.loadFromFile(ocean_grid_file)
    ocean_band = ocean_layer.getBand(0)
    bound = ocean_layer.getBound()  # 获取边界信息

    # 创建结果GridTIFF图层 - 16个波段存储8个分潮的H和G
    band_count = 16
    data_type = DGGS_Float64
    result_layer = DGGSFieldLayerDouble(band_count, data_type, element_type, grid_type, level, bound)

    # 初始化所有网格单元为NAN
    row_size = rowMaxSize(level)
    col_size = fColMaxSize(level, grid_type)
    
    for base_partition in range(1, 11):
        for row in range(row_size):
            for col in range(col_size):
                for band_id in range(band_count):
                    code = DGGSElementCode()
                    code.level = level
                    code.basePartition = base_partition
                    code.elementType = element_type
                    code.row = row
                    code.col = col
                    result_layer.setValue(code, NAN, band_id)

    # 插值数据收集
    if is_interpolate == True and interpolation_count > 0:
        hg_data = {}

    # 收集所有海洋格网ID
    ocean_code_id_list = []
    
    # 遍历所有海洋格网单元
    for base_partition in range(1, 11):
        print(f"处理基础菱形: {base_partition}")
        
        for row in range(row_size):
            for col in range(col_size):
                code = DGGSElementCode()
                code.elementType = element_type
                code.level = level
                code.basePartition = base_partition
                code.row = row
                code.col = col
                
                # 检查是否为海洋格网
                try:
                    ocean_value = ocean_band.getData(code)
                    # 根据您的海洋格网生成代码，海洋格网值为1，陆地为noData
                    if ocean_value != 1:
                        continue
                except:
                    continue
                
                code_id = code2Id(code, grid_type)
                ocean_code_id_list.append(code_id)
                
                if code_id in unique_ids:
                    # 根据id获取所有数据的datatime、ocean_tide
                    code_df = index_df.loc[code_id]
                    # 重置索引
                    code_df_reset = code_df.reset_index()

                    # 选择 datatime 和 ocean_tide 列
                    code_df = code_df_reset[['ocean_tide', 'datatime']]
                    code_df['ocean_tide'] = pd.to_numeric(code_df['ocean_tide'], errors='coerce')
                    code_df = code_df.sort_values(by='datatime')
                    tide_list = code_df.values.tolist()

                    # 转换时间格式
                    elevations = []
                    for item in tide_list:
                        ocean_tide = item[0]
                        datatime = item[1]
                        elevations.append([ocean_tide, NC_DATA_START_DATETIME + datetime.timedelta(days=datatime)])

                    # 最小二乘法
                    lse_result = calLse(wave_names, elevations)

                    if lse_result is not None:
                        result_num = result_num + 1
                        min_hours = (elevations[0][1] - tide_start_time).total_seconds() / ONE_HOUR_TO_SECOND
                        max_hours = (elevations[-1][1] - tide_start_time).total_seconds() / ONE_HOUR_TO_SECOND
                        a0, waves = calHG(lse_result, wave_names, (min_hours + max_hours) / 2)

                        # 按顺序设置16个波段的H和G值
                        band_index = 0
                        for wave_name in wave_names:
                            # 设置H值
                            result_layer.setValue(code, waves[wave_name]['H'], band_index)
                            band_index += 1
                            # 设置G值
                            result_layer.setValue(code, waves[wave_name]['G'], band_index)
                            band_index += 1

                        if is_interpolate is True and interpolation_count > 0:
                            properties = {}
                            for wave_name in wave_names:
                                properties[wave_name + "_H"] = waves[wave_name]['H']
                                properties[wave_name + "_G"] = waves[wave_name]['G']
                            hg_data[code_id] = properties

    # 插值处理
    if is_interpolate is True and interpolation_count > 0:
        inter_hg_data = hg_data
        for i in range(interpolation_count):
            print(f"进行第 {i+1} 次插值...")
            inter_hg_data = interHg(ocean_code_id_list, inter_hg_data, level, grid_type)

        # 将插值结果设置到GridTIFF中
        print("将插值结果写入GridTIFF...")
        for code_id, properties in inter_hg_data.items():
            if code_id not in hg_data:  # 只处理插值得到的格网
                code = DGGSElementCode()
                code.elementType = element_type
                code.level = level
                id2Code(code_id, code, grid_type)
                
                # 按顺序设置16个波段的H和G值
                band_index = 0
                for wave_name in wave_names:
                    # 设置H值
                    result_layer.setValue(code, properties[wave_name + "_H"], band_index)
                    band_index += 1
                    # 设置G值  
                    result_layer.setValue(code, properties[wave_name + "_G"], band_index)
                    band_index += 1

    DestroyKernelObj(kernel)

    # 确保输出目录存在
    output_dir = os.path.dirname(out_file)
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)

    # 保存结果为GridTIFF格式
    result_layer.saveToFile(out_file)
    print(f"计算完成，结果已保存到: {out_file}")
    print("其中共有", result_num, "份可用结果")

    # 释放内存
    del index_df
    del unique_ids


def findInterRange(code, code_ids, grid_type):
    """
    获取插值所需范围

    :param code:        当前格网ID
    :param code_ids:    所有离散格网的ID
    :param grid_type:   离散格网类型

    """
    result = []
    queue = deque()

    kernel = CreateKernelObj(grid_type)
    queue.append(code)
    
    # 插值需要至少4个相邻有效格元
    current = code

    # 找第一圈
    fir_query_result = QueryResult()
    kernel.query(DGGSTopoType_RelativeCells, current, fir_query_result)
    query_result_size = fir_query_result.getSize()
    for i in range(query_result_size):
        temp_code = fir_query_result.getGridCodeByIndex(i)
        query_result_one_id = code2Id(temp_code, grid_type)
        if query_result_one_id in code_ids and query_result_one_id not in result:
            result.append(query_result_one_id)

    # 找第二圈
    if len(result) < 4:
        for i in range(query_result_size):
            current = fir_query_result.getGridCodeByIndex(i)
            sec_query_result = QueryResult()
            kernel.query(DGGSTopoType_RelativeCells, current, sec_query_result)
            sec_query_result_size = sec_query_result.getSize()
            for i in range(sec_query_result_size):
                temp_code = sec_query_result.getGridCodeByIndex(i)
                query_result_one_id = code2Id(temp_code, grid_type)
                if query_result_one_id in code_ids and query_result_one_id not in result:
                    result.append(query_result_one_id)

    DestroyKernelObj(kernel)

    return result


def interHg(ocean_code_id_list, hg_data, level, grid_type):
    """
    对振幅和迟角插值

    :param ocean_code_id_file:  海洋格网文件
    :param hg_data:             计算出来的hg文件
    :param grid_type:           离散格网类型
    """
    no_inter_num = 0
    inter_num = 0

    output = {}
    for code_id in ocean_code_id_list:
        if code_id in hg_data:
            no_inter_num += 1
            continue
        else:
            inter_num += 1
            code = DGGSElementCode()
            code.elementType = DGGSElementType_Cell
            code.level = level
            id2Code(code_id, code, grid_type)

            inter_range = findInterRange(code, set(hg_data.keys()), grid_type)
            properties = {}

            # 初始化所有变量为 0
            K2_H = K2_G = N2_H = N2_G = S2_H = S2_G = M2_H = M2_G = 0
            Q1_H = Q1_G = O1_H = O1_G = K1_H = K1_G = P1_H = P1_G = 0
            for range_code_id_one in inter_range:
                K2_H += hg_data[range_code_id_one]['K2_H']  # K2_H 对应元组第0个元素
                K2_G += hg_data[range_code_id_one]['K2_G']  # K2_G 对应元组第1个元素
                N2_H += hg_data[range_code_id_one]['N2_H']  # N2_H 对应元组第2个元素
                N2_G += hg_data[range_code_id_one]['N2_G']  # N2_G 对应元组第3个元素
                S2_H += hg_data[range_code_id_one]['S2_H']  # S2_H 对应元组第4个元素
                S2_G += hg_data[range_code_id_one]['S2_G']  # S2_G 对应元组第5个元素
                M2_H += hg_data[range_code_id_one]['M2_H']  # M2_H 对应元组第6个元素
                M2_G += hg_data[range_code_id_one]['M2_G']  # M2_G 对应元组第7个元素
                Q1_H += hg_data[range_code_id_one]['Q1_H']  # Q1_H 对应元组第8个元素
                Q1_G += hg_data[range_code_id_one]['Q1_G']  # Q1_G 对应元组第9个元素
                O1_H += hg_data[range_code_id_one]['O1_H']  # O1_H 对应元组第10个元素
                O1_G += hg_data[range_code_id_one]['O1_G']  # O1_G 对应元组第11个元素
                K1_H += hg_data[range_code_id_one]['K1_H']  # K1_H 对应元组第12个元素
                K1_G += hg_data[range_code_id_one]['K1_G']  # K1_G 对应元组第13个元素
                P1_H += hg_data[range_code_id_one]['P1_H']  # P1_H 对应元组第14个元素
                P1_G += hg_data[range_code_id_one]['P1_G']  # P1_G 对应元组第15个元素

            range_length = len(inter_range)
            if range_length < 4:
                print(code_id)
                continue

            properties["K2_H"] = K2_H / range_length
            properties["K2_G"] = float(np.arctan2(
                np.sum(np.sin(np.deg2rad(np.array(K2_G)))),
                np.sum(np.cos(np.deg2rad(np.array(K2_G)))),
            ))
            properties["N2_H"] = N2_H / range_length
            properties["N2_G"] = float(np.arctan2(
                np.sum(np.sin(np.deg2rad(np.array(N2_G)))),
                np.sum(np.cos(np.deg2rad(np.array(N2_G)))),
            ))
            properties["S2_H"] = S2_H / range_length
            properties["S2_G"] = float(np.arctan2(
                np.sum(np.sin(np.deg2rad(np.array(S2_G)))),
                np.sum(np.cos(np.deg2rad(np.array(S2_G)))),
            ))
            properties["M2_H"] = M2_H / range_length
            properties["M2_G"] = float(np.arctan2(
                np.sum(np.sin(np.deg2rad(np.array(M2_G)))),
                np.sum(np.cos(np.deg2rad(np.array(M2_G)))),
            ))
            properties["Q1_H"] = Q1_H / range_length
            properties["Q1_G"] = float(np.arctan2(
                np.sum(np.sin(np.deg2rad(np.array(Q1_G)))),
                np.sum(np.cos(np.deg2rad(np.array(Q1_G)))),
            ))
            properties["O1_H"] = O1_H / range_length
            properties["O1_G"] = float(np.arctan2(
                np.sum(np.sin(np.deg2rad(np.array(O1_G)))),
                np.sum(np.cos(np.deg2rad(np.array(O1_G)))),
            ))
            properties["K1_H"] = K1_H / range_length
            properties["K1_G"] = float(np.arctan2(
                np.sum(np.sin(np.deg2rad(np.array(K1_G)))),
                np.sum(np.cos(np.deg2rad(np.array(K1_G)))),
            ))
            properties["P1_H"] = P1_H / range_length
            properties["P1_G"] = float(np.arctan2(
                np.sum(np.sin(np.deg2rad(np.array(P1_G)))),
                np.sum(np.cos(np.deg2rad(np.array(P1_G)))),
            ))
            output[code_id] = properties

    return output


def calElevation(waves_hg, wave_names, date_time):
    '''
    根据调和常数HG计算平均潮位

    :param waves_hg:    分潮的调和常数
    :param wave_names:  分潮名称
    :param date_time:   计算分潮的时间
    :return:            潮位
    '''
    start_hours = (date_time - tide_start_time).total_seconds() / ONE_HOUR_TO_SECOND  # 当天UTC 0时刻

    # 计算1天24小时的潮位，然后求取平均值
    day_result = np.zeros((HOURS_PER_DAY,), dtype=np.float32)
    for i in range(HOURS_PER_DAY):
        hours = start_hours + i
        wave_fu_dict = cal_f_u(hours)
        wave_v_dict = cal_v(hours)

        result = []
        for wave in wave_names:
            wave_name = str(wave)
            wave_h = waves_hg[wave_name]['H']
            wave_g = waves_hg[wave_name]['G']

            wave_fcosu = wave_fu_dict[wave_name]['fcosu']
            wave_fsinu = wave_fu_dict[wave_name]['fsinu']

            wave_f = math.sqrt(wave_fcosu * wave_fcosu + wave_fsinu * wave_fsinu)
            angle = math.degrees(math.atan(wave_fsinu / wave_fcosu))

            if wave_fsinu >= 0 >= wave_fcosu:
                wave_u = 180 - abs(angle)
            elif wave_fsinu <= 0 and wave_fcosu <= 0:
                wave_u = abs(angle) + 180
            elif wave_fsinu <= 0 <= wave_fcosu:
                wave_u = 360 - abs(angle)
            else:
                wave_u = angle

            wave_v = wave_v_dict[wave_name]
            wave_elev = wave_f * wave_h * math.cos(math.radians(wave_v + wave_u - wave_g))

            result.append(wave_elev)

        result = np.array(result)
        day_result[i] = np.sum(result)

    return float(np.mean(day_result))


def main():

    # 计算调和常数hg
    saveHgFile(
        level, grid_type, element_type,
        input_elev_path, temp_output_path, ocean_grid_path, hg_output_path,
        is_interpolate=True, interpolation_count=4
    )

if __name__ =='__main__':
    start_time = time.time()
    main()
    end_time = time.time()
    execution_time = end_time - start_time
    print(f"程序执行时间: {execution_time} 秒")
